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The multiferroic compound CuO exhibits low temperature magnetic properties similar to an¬ 
tiferromagnetic iron oxides, while the electronic properties have much more in common with the 
high Tc cuprate superconductors. This suggests novel possibilities for the ultrafast optical excita¬ 
tion of magnetism. On the basis on atomistic spin dynamics simulations, we study the effect of 
phonon-assisted multimagnon absorption and photodoping on the spin dynamics in the vicinity of 
the first-order phase transition from collinear to spin-spiral magnetic order. Similar as in recent 
experiments, we find that for both excitations the phase transition can proceed on the picosecond 
timescale. Interestingly, however, these excitation mechanisms display very distinct dynamics. Fol¬ 
lowing photodoping, the spin system first cools down on sub-ps timescales, which we explain as an 
ultrafast magnetocaloric effect. Opposed to this, following phonon-assisted multimagnon excitation 
the spin systems rapidly heats up and subsequently evolves to the noncollinear phase even under 
the influence of isotropic exchange interactions alone. 


I. INTRODUCTION 

The exploration of the speed limits of phase transi¬ 
tions between macroscopically ordered phases of con¬ 
densed matter systems has evolved into a central branch 
of physics, including the study of magnetic phase transi¬ 
tions [1-3], metal-insulator transitions [4-7] and light- 
induced superconductivity [8, 9]. Moreover, recently, 
this investigation has even been extended to multiferroic 
systems, which already in the ground state can exhibit 
both electric and magnetic order and the coupling be¬ 
tween them can give rise to highly non-trivial dynamics, 
such as large amplitude spin-cycloid rotation [10], colos¬ 
sal dynamical magnetoelectric effect [11], and optically 
driven ultrafast transition from collinear antiferromag¬ 
netic (AFM) to spiral states in the multiferroic CuO [12]. 

Gaining microscopic understanding of laser-induced ul¬ 
trafast phase transitions is in general a highly challenging 
computational problem since it requires the modeling of 
the effect of laser absorption, which implies the simula¬ 
tion of coupled time-evolution of the charge, spin and lat¬ 
tice degrees of freedom. Nevertheless, important insights 
can be obtained by breaking this problem into smaller 
parts. For example, laser-induced demagnetization in fer¬ 
romagnetic (FM) metals was successfully modelled by as¬ 
suming that the laser pulse heats the electrons far above 
the Curie temperature on femtosecond (fs) timescales 
[13-18]. Moreover, when different FM metals are coupled 
antiferromagnetically, even laser-heating alone was found 
to be sufficient to cause spin switching [2, 19, 20]. On the 
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other hand, the vast majority of intrinsic antiferromag¬ 
netic materials are insulating, and for excitation below 
the charge transfer gap such ultrafast electronic heating is 
not possible. Nevertheless, ultrafast laser-induced spin- 
reorentation phase transitions have been observed in iron 
oxides, for example, and were explained by the coexis¬ 
tence of two metastable spin-orientation phases enabling 
a spin-inertia driven spin switching mechanism [21, 22] 
and by laser-induced heating of the phonons [23, 24]. 

Here we focus on a different type of spin-reorentation 
transition: the transition from the collinear AFM phase 
to the spin-spiral phase in CuO. Opposed to iron ox¬ 
ides, the electronic structure of CuO has much in com¬ 
mon with the cuprate high-temperature superconductors, 
yielding different possible excitation mechanisms. First, 
in these materials it is well-known that the optical ab¬ 
sorption below the charge transfer gap is dominated by 
phonon-assisted multimagnon (MM) excitation [25, 26]. 
Second, for these Mott-insulating materials it was re¬ 
cently predicted that photo doping (PD) can cause an 
ultrafast modification of the exchange interaction, with 
an effect comparable to chemical doping [27] . The latter 
is particularly interesting since recent investigations show 
that chemical doping dramatically reduces the tempera¬ 
ture at which the antiferromagnetic to spin-spiral phase 
in CuO occurs [28]. Hence, both MM excitation and PD 
suggest distinct and rather direct ways to induce this 
phase transition, which so far have not been explored. 

In this article we investigate systematically phonon- 
assisted MM excitation and PD as mechanisms for the 
laser-induced phase transition in CuO on the basis of 
atomistic spin dynamics (ASD) simulations. We demon¬ 
strate that both excitation mechanisms can induce the 
phase transition on picosecond (ps) timescales. Intrigu- 
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ingly, we find that different excitation mechanisms nev¬ 
ertheless induce qualitatively very different dynamics in 
which the spin system either cools down or heats up on 
sub ps timescales. In addition, we find that the fastest 
phase transition occurs for phonon-assisted MM excita¬ 
tion and demonstrate that only in this case the phase 
transition proceeds even in the absence of anisotropy and 
dissipation. This suggest that exchange interactions be¬ 
tween the spins alone are sufficient to drive the system 
across the transition, which explains why it can be faster 
than the dynamics caused by PD. 

The paper is organized as follows: The properties of 
CuO are discussed in Sec. II. The model for CuO and 
the semiclassical equation of motion method are intro¬ 
duced in Sec. III. The different excitation mechanisms 
are treated in Sec. IV and we discuss how they are mod¬ 
elled within the ASD simulations. The results from such 
simulations are presented in Sec. V and compared with 
the results reported from ultrafast time resolved pump- 
probe measurements [12]. A discussion and outlook is 
given in Sec. VI. 


II. MAGNETIC PROPERTIES OF CUO 

CuO is a Mott insulator and spin S' = 1/2 compound 
with a strongly dominant exchange interaction along the 
crystallographic [101] axis. At high temperatures it is 
effectively a one-dimensional quantum spin chain with 
fermionic spinon excitations [30]. Below a weakly first 
order transition at Tjv 2 = 230 K, a three dimensional 
magnetic ordering (AF2) sets in, which over a 17 K tem¬ 
perature span down to Tati = 213 K, is an incommensu¬ 
rate spin ordering with nearest neighbour spins at nearly 
90 degrees to each other [Fig. 1(b)]. This phase has been 
identified as multiferroic by Kimura and collaborators 
[31]. 

Below the first order phase transition at Tati, a com¬ 
mensurate collinear ordering (AFl) is found [Fig. 1(a)]. 
The ordering of the AFl and AF2 phases and the overall 
temperature dependence of the phase diagram [Fig. 1(c)] 
can be understood in terms of a competition between 
the frustration among the Heisenberg couplings, which 
promotes the AF2 configuration, with the biquadratic 
exchange couplings, which suppresses the non-collinear 
ordering for temperatures up to Tati [28, 32]. 

The unusually high temperature of the AF2 multifer¬ 
roic phase has been a strong driving force for experi¬ 
mental and theoretical activity [29, 33-40] on cupric ox¬ 
ide, CuO. In this phase the canted magnetism drives im¬ 
proper ferroelectricity[28, 29, 33-36]. Whereas the un¬ 
derstanding of the equilibrium phase diagram has pro¬ 
gressed, it is unclear how the phases evolve into each 
other upon laser excitation. 

Using a free electron laser as source for the probe beam, 
time-resolved diffraction experiments revealed an ultra¬ 
fast optically driven phase transition in CuO from the 
AFl phase to the AF2 phase [12]. The diffraction signal 
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AF1 phase: 0 < T < 213 K AF2 phase: 213 < T < 230 K 



FIG. 1. The spin configurations of the (a) AFl and (b) AF2 
phases in CuO. Blue and grey spheres are Cu atoms on differ¬ 
ent (010) planes [28, 29]. The red arrows show the Cu spin po¬ 
larization. Entering the AF2 phase, the qAFi = (0.5, 0, —0.5) 
ordering vector is modulated to the incommensurate order¬ 
ing vector qAF 2 = (0.506, 0, —0.483). (c) The static structure 
factor S'(q) for qAFi and qAF 2 , and the total energy E and 
heat capacity Cv per spin, obtained in classical Monte Carlo 
simulations for a 36a x 4b x 36c supercell using the parame¬ 
ters of Ref. [28]. The structure factor is normalised to unity 
for the T = 0 K AFl spin configuration. 


gave immediate insight into the evolution in time of the 
magnetic ordering upon excitation. Measurements were 
performed at T = 207 K, slightly below Tati. The pro¬ 
cess was fluency-dependent up to a saturation fluency for 
which a minimum transition time of 400 fs was reported. 
In the next section we discuss the model and methods we 
use to simulate the dynamics between these two phases. 


III. MODEL AND METHODS 

The phase competition between the AFl and AF2 
phases in CuO can be captured by a classical spin 
Hamiltonian which includes the competition between bi¬ 
linear Heisenberg exchange and biquadratic exchange 
[28, 32, 36]. 

In recent theoretical studies, density functional theory 
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calculations have been pursued to parametrise these ex¬ 
change interactions from the total energies for different 
spin conhgurations [28, 33, 35, 41]. For the present in¬ 
vestigations the Hamiltonian 

Jtyi = ^xch + + ^ni U j 

= ^ E ■ S, + i E i E 

is used, where the first term is Heisenberg exchange 
(exch), the second term is biquadratic exchange (bq) and 
the third term is anisotropic exchange (ani) and is a 
classical spin of magnitude |Si| = Si = 1/2. 

With the set of parameters reported in Ref. [28] classi¬ 
cal Monte Carlo simulations of Eq. 1 produce a phase dia¬ 
gram with a first order phase transition at Tni ~ 165 K, 
and a weakly first order phase transition to the para¬ 
magnetic phase at Tn 2 ~ 179 K. In Figure 1(c) is 
shown the structure factor S'(q) for the ordering vectors 
Qafi = (0.5,0,—0.5) and qAF2 = (0.528,0,-0.472) in a 
temperature interval around the phase transition at Tni . 
To go from the commensurate AFl phase to the incom¬ 
mensurate AF2 phase latent heat needs to be supplied 
represented by the sharp peak in the specific heat. Ex¬ 
cept for the underestimation of the Neel temperatures, 
the phase diagram is in excellent agreement with experi¬ 
ment for the nature of the phases, wave-vector and nature 
of the order parameters, etc [28]. 

Given the good agreement for the equilibrium phase di¬ 
agram, we will use ASD simulations based on the Hamil¬ 
tonian Eq. (1) as minimal approach for studying the 
nonequilibrium spin dynamics across the AF1-AF2 phase 
transition. Within ASD, the dynamics of a spin sys¬ 
tem is modelled seiniclassically by the stochastic Landau- 
Lifshitz equation (SLL) 

= -S, X [B, + Bf (t)] - X {S, X [B, -h Bf (t)]}, 

( 2 ) 

which describe the motion of the classical spins in an 
effective magnetic field Bi, calculated from Bi = — §^. 
Bf (t) is a stochastic magnetic field with a Gaussian dis¬ 
tribution. By the fluctuation-dissipation theorem, the 
magnitude of B[* (t) is related to the dimensionless damp¬ 
ing parameter a. 

We note that the dynamics we investigate are at tem¬ 
peratures high enough to allow for the assumption of a 
classical renormalised regime [42] which renders Eq. (2) 
applicable even for a spin 1/2 system. 

To characterise the dynamics across the phase tran¬ 
sition we focus on two observables. First, we com¬ 
pute the equal time correlation function S'(q, t) = 
where C{ri - rj,t) = [S,(t) • Sj(t)] 
is the spin-spin correlation function for spins displaced 
by the vector — Vj. [...] indicate averaging over the 
cell and over different realisations of the initial spin ex¬ 
citations and the heat bath. S{q^,t) shows two peaks at 
q = qAFi and q = qAF 2 which we use to monitor the 
evolution from one phase to the other. 


Second, we monitor the evolution of the temperature 
of the spin system. To this end, we sample the energy Ei 
defined as the difference 

E, = E- Ef, (3) 

between the instantaneous total energy E and the total 
energy Ef when spin is relaxed keeping all other spins 
fixed. The energy Ef is minimised, for fixed orientation 
of the other spins with a gradient descent algo¬ 

rithm. The temperature T is then obtained by fitting 
the histogram of {Ei} to the Boltzmann distribution 

P[Ei\=Pi,e-f^^\ (4) 

where ks is the Boltzmann constant and /3 = l/iksT) 
is the inverse temperature. The method of fitting the 
histograms of {Ei} to a Boltzmann distribution assumes 
that the spins are close to the local equilibrium state 
and is therefore only applicable for moderate tempera¬ 
tures. We have investigated the range of applicability 
in more detail in Appendix A, including a quantitative 
comparison with an alternative approach for computing 
the spin temperature [43] , and found our approach to be 
sufficiently accurate to capture the temperature regime 
studied below. 


IV. EXCITATION MECHANISMS 

In general, the light-matter interaction in the optical 
regime relies on the electric dipole approximation. This 
makes the direct excitation of the magnetic degrees of 
freedom in materials with inversion symmetry weak or 
even forbidden. In addition, different from magnetic met¬ 
als, in which the rapid heating of the electrons turns out 
to be very efficient for driving (sub) picosecond phase 
transitions [19, 20], CuO is insulating and photon absorp¬ 
tion does not seem to be relevant on first sight. Below, 
however, we discuss three distinct mechanisms which in¬ 
directly couple light to the magnetic degrees of freedom. 
All mechanisms rely, on one way or the other, on modu¬ 
lation of exchange interactions J, either by joint excita¬ 
tion of a phonon (phonon assisted infrared absorption), 
by changing the electronic distribution (photodoping) or 
by modulating J with an electric field (two-magnon Ra¬ 
man scattering). We then estimate their efficiency based 
on literature values and explain how we effectively in¬ 
clude these excitation mechanisms within atomistic spin 
dynamics simulations. 

A. Infrared absorption 

In materials with inversion symmetry, such as CuO in 
the AFl phase, the excitation of magnons by absorption 
of infrared photons is in principle electric dipole forbid¬ 
den. However, the inversion symmetry can be broken by 
the excitation of a phonon which enables simultaneous 
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excitation of an even number of magnons [25, 26]. This 
mechanism dominates for infrared absorption in CuO lay¬ 
ers and similar magnetically ordered insulators and is de¬ 
scribed by the Hamiltonian: 

H = 

P ph E -b Hph (5) 

i,5 


Here E describes the electric field of light, Ui the phonon 
displacements and J{E, Ui+ 5 / 2 ) the exchange interaction 
modified by the electric field and the phonon displace¬ 
ments. Further S indicates nearest neighbor sites with 
respect to site i, while Pph = J2i,s ^oUi+s /2 and Hph 
are the phonon dipole moment and phonon Hamiltonian, 
respectively. 

Expanding the Hamiltonian to linear order in the elec¬ 
tric field and the ionic displacement one obtains the 
dipole moment for Cartesian component 


dJ 
dJ 


Pu = P, 

E 


SiSi^S 


( 6 ) 


i.Sa 


dEpdui_^_s/2, 


/2,i^ ^ i+5 


Neglecting the small ferroelectric effects in the multifer- 
roic phase the second term is zero in the presence of in¬ 
version symmetry. The third term is the one responsible 
for the phonon-assisted absorption of magnons. 

We observe that the phonon breaks the inversion sym¬ 
metry as required for non-vanishing electric-dipole tran¬ 
sitions. We note that opposed to the AFl phase, the 
AF2 phase exhibits a ferroelectric dipole moment which 
in principle facilitates a direct coupling of the electric 
field component of the laser pulse to the magnetic degrees 
of freedom through the second term of Eq. (6). However, 
such coupling is only effective once the system is already 
in the AE2 phase, whereas we want to address how fast 
the reorientation of the spin system from the collinear 
AEl phase to the incommensurate spin spiral AF2 phase 
can occur. In any case, we expect spontaneous inversion 
symmetry effects to be much smaller since the involved 
dipoles are minute. Indeed, they are beyond the sensi¬ 
bility of X-ray experiments, for example, and are only 
detected through the macroscopic polarization. 

Clearly, it would be very interesting to study the com¬ 
plete dynamics of the spin and lattice degrees of free¬ 
dom that follow from Eq. (5). This, however, is very 
challenging since it requires at least to go beyond the 
semi-classical description for the spin dynamics. In ad¬ 
dition, augmenting classical phonon dynamics and their 
coupling to ASD would only slightly change the spin dy¬ 
namics in the regard that displacement (spin) correla¬ 
tions renormalize the magnon (phonon) frequencies, but 
this spin-phonon coupling is not strong enough to change 
the order of time-scales of the transversal spin reorienta¬ 
tion. Therefore, here we focus exclusively on the classical 
spin dynamics following phonon-assisted MM excitation 
and omit solving directly the time evolution during the 
excitation of magnons and phonons. 


To effectively take this excitation into account we as¬ 
sume each absorbed photon creates a phonon plus mul¬ 
tiple two-magnon excitations. We neglect the first and 
simulate the second by an exchange process. Namely, 
we interchange a fraction of nearest neighbor spins along 
the [10l]-direction, which has the dominant exchange 
coupling in CuO [28, 30]. To estimate the fraction of 
spin pairs that are interchanged, we use that in middle- 
IR absorption measurements with photon energies be¬ 
low the charge transfer gap, the absorption coefficient of 
CuO was determined to be of the order a ^ 200 cm“^ 
(Ref. [44]), with corresponding skin penetration depth 
d ^ 0.005 cm. Defining the average number of absorbed 
photons per unit volume A^phot, the number p of absorbed 
photons per Cu atom can be computed through the ex¬ 
pression 


P — A^photEcu 


Apump 

Aphot 



( 7 ) 


with p taking values 4.6-10“^ (1.8-10“^) for laser fluency 
7 mJ/cm^ (28 mJ/cm^). Each absorbed photon cause 
multiple phonon assisted magnons pairs [25, 26, 45]. 
Therefore, in the simulations that are presented in Sec. V, 
a fraction in the range [a; = 0.005,... 0.040] of randomly 
chosen pairs of spins are interchanged. This yields a 
small but finite probability also for four magnon scatter¬ 
ing which take place when two pairs of exchanged spins 
are in the immediate vicinity of each other. 


B. Photo doping 

Photon absorption at frequencies comparable to the 
gap can cause PD, in which the electronic distribution 
is changed, which subsequently can modify the exchange 
couplings and spin ordering [3, 27]. For Mott-insulators 
like CuO, PD comprises the creation of holon-doublon 
pairs (h-d pairs) and their relaxation within the Hub¬ 
bard bands causes a rapid reduction of the ordered spin 
moment (S^) and the exchange parameter within a few 
tens of the hopping time [27]. Since electron hopping 
occurs on the fs time scale, PD has the potential of be¬ 
ing fast enough for launching spin dynamics that lead to 
the phase transition on (sub) picosecond time scales. In 
a classical picture, a simplified way to account for the 
presence of h-d pairs is by introducing vacancies. Alter¬ 
natively, one can introduce reduced local moments [sim¬ 
ilarly to a reduction of S'(qAFi,i)], in conjunction with 
reduced exchange parameters. For small concentrations 
X we have checked these approaches produce qualitatively 
very similar results. In the simulation results presented 
below, we have therefore restricted ourselves to the in¬ 
troduction of vacancies. Note that, we could also let the 
next nearest neighbor (NNN) couplings of the vacancy 
be changed, in line with the second order Henley effect 
in Ref. [28] , but for simplicity we limited ourselves to the 
introduction of vacancies alone. In addition, we assume 
that the vacancies persist over time scales longer than 
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the ASD simulation time (t > 5 ps), in order to isolate 
purely the effect of PD alone. For the concentration of 
vacancies, values in the range [x = 0.001,.. .0.050] have 
been used (c.f. the case of phonon assisted MM excita¬ 
tion). We concentrate on short times and assume a long 
life time of the holon-doublon pairs. The later decay of 
the holon-doublon pairs through MM decay [46-48] con¬ 
stitutes another channel for driving the spin dynamics, 
which we neglected for simplicity. 

C. Two-magnon Raman scattering 

For comparison to the above mentioned absorption 
mechanisms, we also discuss briefly the two-magnon Ra¬ 
man scattering. Here an incoming photon scatters inelas- 
tically in an interaction where two magnons with oppo¬ 
site momenta are created [49]. For high Tf. cuprate su¬ 
perconductors, the assignment of the spectral weight to 
two-magnon scattering has been used to determine the 
values of the leading exchange interactions in the CUO 4 
planes [50]. For CuO, two-magnon scattering was dis¬ 
cussed as contributions to the broad peak features around 
1100 and 2000 cm“^ [51]. The fraction of spins which are 
excited in CuO through two-magnon Raman scattering 
during a pump pulse can be estimated from the cross- 
sections known from spontaneous Raman scattering on 
chemically related copper oxides [52]. For YBa 2 Cu 306 .o 
it was reported [52] a two-magnon Raman cross section 
<S' 2 M = 1-2 ■ 10 “^ cm“^sr“^ per unit volume of sample. In 
the pump probe experiment of [ 12 ] the pump laser wave 
length A = 800 nm corresponds to a photon energy of 
Aphot = 1-55 eV. The intensity of the pump laser was 
such that over the 40 fs FWHM pulse lengths, a fluency 
of ifpump = 7, 14 or 28 mJ/cm^ were incident on the 
sample. The 4 atom primitive cell of CuO has volume 
Vpc = 0.081 nm^. The total number of photons per unit 
area impinging the surface of the sample during one pulse 
is given by A^phot = Apump/Aphot. Multiplying by S 2 M 
and by the volume per Cu, and integrating over 47 r sr we 
obtain the total number of spin flips created per Cu. 

p = 47rS'2MA^photVcu = 47rS'2M^^^Vcu (8) 

-^phot 

For a fluency ifpump = 7 mJ/cm^ a fraction of two- 
magnon spin flips p ~ 10 “® is obtained, a value too small 
to allow inelastic scattering to be of relevance for driving 
the spin system from the AFl to the AF2 phase. There¬ 
fore, below we restrict ourselves to the spin dynamics 
that follows from the absorption of light. 

V. SIMULATION RESULTS 

In this section we present simulations results in which 
spin system is excited with a varying a fraction of inter¬ 
changed spin-pairs and photo-induced vacancies, simulat¬ 
ing phonon-assisted MM excitation and PD, respectively. 


In addition, for comparison we simulate the response to 
sudden changes of the temperature, mimicking phonon¬ 
heating. For each case, we have followed the evolution 
in time of iS'(q, t). All simulations have been performed 
with the UppASD software [53] for a 36a x 4b x 36c 
supercell with periodic boundary conditions. For this 
cell size S'(q, t) has peaks for qAFi = (0.5,0,0.5) and 
<S'qAF 2 = (0.528,0, —0.472) and we use the 5'(q, t) values 
for these wave vectors to follow the time evolution of the 
magnetic ordering. We approximate the excitations to 
be instantaneous, i.e at time t = 0 we swap a fraction x 
of pairs of spins or include a fraction x of vacancies in 
the system. This is a simplification of the actual opti¬ 
cal excitation process which takes place for an incident 
laser pulse that is few tenths of femtoseconds long and 
the photon-matter interaction developing over multiples 
of electronic hopping times. 

In order to clarify the microscopic origin of the ob¬ 
served dynamics, below we will present simulations with 
and without coupling to the heat bath. The latter gives 
undamped (a = 0) evolution. Under the isotropic Hamil¬ 
tonian + both angular momentum and 

total energy are constants of motion, but energy redistri¬ 
bution between different degrees of freedom can still give 
rise to non-trivial dynamics. At finite damping, a > 0, 
energy can dissipate to the environment, enabling the 
system to relax back to the initial state. 

We also note that, despite CuO being an antiferro- 
magnet, at finite temperature there is anyhow a small 
but finite net magnetic moment, given the circumstance 
that the thermally fluctuating spins do not fully com¬ 
pensate each other to give identically zero total magnetic 
moment. In equilibrium at T = 160 K, the magnetic mo¬ 
ment in Langevin dynamics a = 0.01 at times t < —10 ps 
fluctuates with an amplitude ~ 10 “^^b- Using the semi- 
implicit solver SIB [54] for the Landau-Lifshitz equation, 
the total magnetic moment and its vector components 
are conserved to ±10“^5 in a = 0 evolution over the 
simulation interval —10<t<0ps. The total energy 
of the spin system is conserved to ±10“^ meV over the 
same simulation time. 

The protocol for the simulations has been the fol¬ 
lowing: Starting from the T = 0 K AFl configuration 
the spins are relaxed to be in equilibrium with a fi¬ 
nite temperature. The system then evolves for times 
— 10 < t < 0 ps, to check that spontaneous AF1-AF2 
transitions do not occur, before it is excited at t = 0 
and subsequently evolved to extract the spin dynamics 
caused by photo-excitation. 


A. Phonon-assisted multimagnon excitation 

In this section results are reported from simulations 
of phonon assisted MM excitation, which dominates for 
photon energies below the charge transfer gap. First, 
we focus on simulations at low temperature T = 10 K, 
to show that this causes a nonequilibrium distribution 



6 


of the spins which are elucidated by histograms of {Ei} 
energies. Second, we show simulation results for the ini¬ 
tial temperature T = 160 K, from which MM excitation 
turns out to be sufficient to drive the system in the AF2 
phase. 

In Figure 2 is shown the a = 0 (full lines) and a = 0.01 
(dashed lines) evolution in time of the energy histogram 
after an excitation of x = 0.02 (the concentration of pairs 
of spins exchanged) of two-magnon excitations of a spin 
system initially in thermal equilibrium at T = 10 K. 
Upon excitation, energy is injected to the spin system 
which is brought to a non-equilibrium state. The black 
symbols display the energy distribution before pump 
when the spin system is in equilibrium. The red lines 
for f = O’*' fs (after pump) have distinct peaks (marked 
out with blue arrows) indicating the energy and the in¬ 
tensities of two- and four-magnon excitations. The in¬ 
sets show the evolution of the temperature of the sys¬ 
tem, with the first few data points displaying least mean 
square values of the temperature in a non-equilibrium 
regime in which temperature is not well defined. After 
t = 1.0 ps the system has equilibrated to a new tem¬ 
perature of T = 50 K. We note that this relaxation to 
a Boltzmann distribution took place in microcanonical 
evolution. For an isolated spin system this kind of equi¬ 
libration can occur due to the non-linearities inherit in 
the Landau-Lifshitz equation, but it is not possible for 
e.g. for an integrable system like an isolated system of 
coupled harmonic oscillators. 

For the case of evolution a = 0.01 with finite damping 
(dashed lines), the excitation energy added to the system 
dissipates as the system relax to be in equilibrium with 
the heat bath. For times up to t = 150 fs, the evolution 
of the energy distribution is very similar to the system 
with undamped dynamics. At t = 1.0 ps, the equilibra¬ 
tion is essentially complete and the system has reached 
a temperature of T = 10.0 K. 

We emphasize that at low temperature distinct peaks 
corresponding to two- and four-magnon excitation can be 
seen, indicated by arrows in Figure 2. Below we study 
MM excitation at much higher temperatures. In that 
case, there is already a large fraction of spins that have 
energies that are of the same magnitude as the energy of 
a MM excitation event and therefore such sharp features 
are absent. From this we understand that the average 
energy given to the spins by MM excitation is in itself 
temperature dependent. It will take a higher value at 
low temperatures where the spins are close to perfectly 
collinear, but a reduced value at intermediate and high 
temperatures where the effective fields are reduced by 
thermal fluctuations. 

The X = 0.02 excitation did not bring enough energy 
to push the r = 10 K system into the AF2 regime, 
which over few ps relaxed to 50 (10) K for undamped 
(damped) dynamics. The situation changes if the exci¬ 
tation takes place closer to Uni- In Figure 3 is shown 
the results of a simulation where the spin system, ini¬ 
tially in thermal equilibrium at T = 160 K is excited 



FIG. 2. Histogram of energies {Ai} after exciting a system in 
equilibrium at T=10 K (black circles) with a concentration 
X = 0.02 of two-magnons flips. The blue arrows indicate 
the two peaks in the red curve that are corresponding to 2 
and 4 magnon excitations respectively. The inset shows the 
temperature evolution of the spin system. Full (dashed) lines 
show undamped a = 0 (damped a = 0.01) dynamics. 


with an excitation intensity of x = 0.02 and evolved 
in undamped a = 0 (full lines) and damped a = 0.01 
(dashed lines) dynamics. For a = 0 the system under¬ 
goes a phase transition from AFl to AF2 and reaches 
new equilibrium values of the order parameters after ap¬ 
proximately 4 ps. An initial decay can be observed for 
and also for S{q 2 ,t), due to the circumstance 
that the sublattice magnetization is reduced abruptly in 
the two-magnon excitation. This can be observed more 
directly in the inset, in which the normalised quantities 
S'(qi, t)/S'(qi, t = 0) and S'(q2, t)/S'(q2, t = 0) are plot¬ 
ted. We note that within the insets time window of 0.25 
ps, the a = 0 and a = 0.01 dynamics depart only slightly 
from each other. As expected the decrease of S(qi,t) is 
monotonous, but also S(q 2 ,t) decay over the first 50 fs. 
For times t > 50 fs, S(q 2 ,t) is increasing as the spin 
system enters the transition into the AF2 phase. 

For a system in equilibrium at T = 160 K, the change 
to the shape of the energy histogram on two-magnon exci¬ 
tation is much less drastic as compared to the T = 10 K 
case shown in Fig. 2. We visualize the deviation from 
the thermal equilibrium distribution /„(t < 0) by plot¬ 
ting the relative histogram /„(t)//„(f < 0), where is 
the bin count of the n:th bin, as shown in Fig. 5 for 
a = 0 [Fig. 5(a)] and a = 0.01 [Fig. 5(b)]. This clearly 
demonstrates that, similar as in the T = 10 K case, MM 
spin flips cause high energy excitations in the spin sys¬ 
tem. The quota /„(t)//„(f < 0) becomes very noisy at 
at energies larger than 35 meV, given the overall low bin 
counts in this energy range. Therefore /„(t)//„(t < 0) is 
shown for the energy interval 0 to 35 meV. 
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time (ps) 


FIG. 3. Dynamics following multimagnon excitation. The 
upper panel shows the trajectory in time of ^(qi, t) (red) and 
S(q 2 ,t) (blue), averaged over 128 replicas, for a = 0 (full 
line) and a = 0.01 (dashed line) dynamics, t < 0 is before 
pump and t > 0 is after pump. The initial temperature is 
T — 160 K, and the initial excitation intensity is a; = 0.02. 
The lower panel shows the evolution of the temperature of 
the system. The inset shows the first 250 fs of evolution for 
the normalised order parameters. The dash-dotted lines are 
for simulations using the isotropic Hamiltonian Jiflso- 


The lower panel of Figure 3 showns the evolution in 
time of the temperature which for the case of undamped 
dynamics initially peaks at T = 170 K. In the subse¬ 
quent evolution the temperature decreases to T = 169 K 
as latent heat is necessary to enter the AF2 phase. As 
can be seen in Fig. 4, the energy injected by the initial 
excitation brings the average spin energy to the value 
which in thermal equilibrium corresponds to a temper¬ 
ature T Rs 169 K. For a = 0.01 the order parameters 
and 5'(q2,t) display a transient reduction and 
enhancement respectively, and the system relaxes back 
to T = 160 K within 4 ps. 

In order to analyse how large influence anisotropic 
terms in the Hamiltonian have on the phase transition, 
we have pursued simulations for the isotropic Hamilto¬ 
nian J^iBo ■ The instantaneous excitation of the spin sys¬ 
tem at t = 0 in form of swaps of neighboring spins, inject 
energy but conserve the angular momentum since the 
swap operator is rotationally invariant. The dash-dotted 
lines in Fig. 3 show the evolution in time of 5'(qi, t) and 
S{q 2 ,t) respectively in simulation of Jtfso- Interestingly, 
we find that the transition between the AFl phase and 
the AF2 phase after the initial injection of energy by MM 
excitation, can proceed solely through redistribution of 
energy between the isotropic degrees of freedom. 

Johnson et al. [12] investigated the dependence of a 
switching delay time defined as the onset time at 
which the ratio S'(q2,t)/5'(qi,t) would start to grow, 
on the pump fluency. Starting from lower fluencies 



FIG. 4. (Golor online) The full line shows the energy of the 
spin system as a function of temperature. The blue and black 
symbols indicate the initial T and the energy after excitation 
with strength x =0.010 or 0.020 respectively (the concentra¬ 
tion of pairs flipped) of a system which before pump is in 
thermal equilibrium at T — 160, 162 or 164 K. 


decreased with increasing fluency, but it was not pos¬ 
sible to push it below = 400 fs. In the present 
model for phonon assisted MM excitation, the analogue 
of pump fluency is the concentration of initial spin pairs 
flipped. Results from simulations with initial concentra¬ 
tion in the range [0.005, 0.040] are presented in Fig. 6 
and we make the following observations; i) In the inter¬ 
val 100 <t< 400 fs, the rate of growth of S{q 2 ,t) is very 
similar in the fluency ranges [0.010, 0.040], with x = 0.020 
being the concentration that for t = 5 ps has given the 
highest value of 5'(q2,t). For fluencies x = 0.030 and 
0.040 the spin system is over pumped and brought closer 
towards the paramagnetic region, ii) The inset displays 
the evolution for the first 250 fs and reveal that the time 
tp™, at which S{q 2 ,t) start to grow after the initial dip, 
takes the value = 50 fs independent of fluency. 

For the excitation x = 0.02 it can be observed in 
the inset of Fig. 3, that the onset time t®™ is not af¬ 
fected by damping a = 0.01. This can be understood 
from the following considerations: With the AFl as ref¬ 
erence spin configuration, the highest precession frequen¬ 
cies occur when summing up the exchange interactions 
to Y.i Jij = 168 meV for /max = T^max = 81.8 THz, 
which gives for a = 0.01 the corresponding frequency of 
damping motion /“ax = a/max = 0.818 THz, or a mini¬ 
mum relaxation time Tmin = l/(27r/“ax) = 194 fs. Thus, 
damping motion will give significant influence only after 
hundreds of fs, but will for a = 0.01 be negligible for the 
first few hundred fs of evolution, constituting a regime of 
exchange interaction driven dynamics [19, 55]. 
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FIG. 5. The relative histograms < 0) of energies {Ei} 

evolving in a) a = 0 and b) a = 0.01 dynamics respectively, 
after exciting a system in equilibrium at T=160 K (black cir¬ 
cles) with a concentration x = 0.02 of two-magnons flips. The 
quota In{t)/In{t < 0) becomes very noisy at at energies larger 
than 35 meV, given the overall low bin counts in this energy 
range. Therefore In{t)/In.{t < 0) is shown only for the en¬ 
ergy interval 0 to 35 meV. The insets show the temperature 
evolution of the spin system. 


B. Photodoping induced dynamics 


The energy of the incident 800 nm laser pump pulse 
used in [12] coincide with the charge transfer (CT) gap 
of CuO. In this regime one can expect that in addition to 
the previous mechanism also PD to be an important exci¬ 
tation mechanism. Photodoping will have two effects. If 
the incident energy is larger than the CT gap the excited 
particle hole pairs will relax by the emission of magnon 
and phonon excitations until they reach the minimum 
possible unbound CT energy or an excitonic bound state. 
This transfer of energy can be modelled as a heating of 
the phononic system or the magnetic system and in our 
simulations corresponds to a rise of the bath tempera¬ 
ture in the damped dynamics. On the other hand, as we 
will show below, there is a magnetocaloric effect due to 
the sudden change in the magnetic Hamiltonian. We will 
analyze this effect in detail neglecting the heating effect, 
which can be minimized by tuning the incoming laser 
to the minimum CT energy. The temperature induced 
dynamics is discussed below. 

As discussed in detail in Ref. [28] , in equilibrium one of 
the main effects of dilution of the spin system is to shift 
the magnetic phase diagram to lower temperatures. Be¬ 
low we present nonequilibrium studies in which we follow 
the dynamics after sudden dilution. The concentration 
of vacancies x is chosen in a similar range as the pairs 
of interchanged spins discussed with results of phonon 
assisted MM excitation. 



FIG. 6. Excitation intensity dependence of dynamics follow¬ 
ing MM excitation. The trajectories in time of S(qi, t) (upper 
panel) and S{q 2 ,t) (lower panel), averaged over 128 replicas, 
and with a = 0. The initial temperature was T — 160 K, and 
the initial excitation intensity was in the range [0.005, 0.040]. 
The inset shows the evolution of the order parameters during 
the first 250 fs. 


In Fig. 7 is shown the evolution of the order param¬ 
eters ^(qi,!) and S{q 2 ,t) following on PD that intro¬ 
duces X = 0.02 vacancies. Full lines show evolution at 
a = 0. The S{qi,t) ordering experience an abrupt step 
as the sublattice moment is reduced, correlated with the 
reduced number of spins in the system. After this initial 
drop, S'(qi,t) levels out to a slightly reduced value and 
S'(q 2 ,t) is, over the time scale of 5 ps, not affected. In¬ 
terestingly, the temperature of the system experiences a 
rapid drop from T = 160 K to 157 K. We explain this 
drop as an ultrafast magnetocaloric effect, where the role 
of the external magnetic field is replaced by the change of 
the effective exchange field on the spins by photo-doping. 

To substantiate our interpretation of a magnetocaloric 
effect we computed the equilibrium entropy per spin as 
a function of T for different vacancy concentrations x, 
S{T,x) = So{x) + /^/rCM(r, x) jr and the magnetic 
heat capacity Cm is computed from the energy fluctu¬ 
ations {{5EY). Here the total entropy is divided by 
the actual number of spins in the sample, not the num¬ 
ber of sites. To obtain an accurate account of the en¬ 
tropy differences, we matched So{x) to 5*0(0) in the high- 
temperature paramagnetic phase at T = 1260 K, while 
To = 1 K. 

The curves shown in Fig. 8 are qualitatively similar to 
what would be obtained for independent spins in a uni¬ 
form magnetic field, with the upper curves corresponding 
to lower magnetic field. The reason is that both the mag¬ 
netic field and the interaction constrain the spin motion, 
effectively reducing the fluctuating phase space and low¬ 
ering the entropy. Taking a reference spin, after a sudden 
quench of the magnetic field (classical magnetocaloric ef¬ 
fect) or after a neighbor is eliminated (present Simula- 
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tion) the spin finds itself in a state which is much more 
ordered than what it should be at the initial tempera¬ 
ture. This corresponds to the adiabatic jump shown in 
Fig. 8 with the dashed line. The vertical arrows show 
the temperatures obtained in spin dynamics simulations 
after internal equilibration in a = 0 dynamics, with the 
error bars stemming from the temporal fluctuations and 
the measurement of the temperature (see Appendix A). 
We can observe that in the spin dynamics simulations the 
temperatures are lowered almost but not quite as much 
as the constant entropy jump implies. 

Similar rapid drops in the temperature are obtained 
with changing the strength of the Heisenberg exchange 
interaction {1 — x)Jz (data not shown). In this 

case we can estimate the effective field seen by a sin¬ 
gle spin simply from zxJzS'^ = gfiBSBgS, with the 
dominant Heisenberg exchange interaction [28], yielding 
BeS = 10 T at X = 0.01 and an efficiency of 0.15 K/T. 
This is comparable but slightly lower than the efficien¬ 
cies reported in [56] and references therein, which is un¬ 
derstandable since our results are obtained deep in the 
ordered phase of CuO. 

Although the system thermalizes even at a = 0, 
S (q 2 , t) is not affected on the time scale of the simu¬ 
lation. Hence, despite the fact that the introduction of 
vacancies reduces T^vi, our simulations do not show the 
evolution to a different phase. We can understand this 
since the shift in Tjvi is smaller than the magnetocaloric 
effect. Indeed, already at x = 0.02 we obtain a reduc¬ 
tion to T Ri 156 K from the magnetocaloric effect, while 



time (ps) 


FIG. 7. Dynamics following on PD. The upper panel shows 
the trajectory in time of S{qi,t), averaged over 128 replicas, 
for Q = 0.00 (full lines) and a = 0.01 (dashed lines). The 
initial temperature was T = 160 K, and the initial excita¬ 
tion intensity of vacancies was x = 0.02. The lower panel 
shows the evolution of the temperature of the system.The in¬ 
set shows the first 250 fs of evolution for the normalised order 
parameters, using the same range of the axes as in the inset 
of Fig. 3. 



FIG. 8. The entropy S{T,x) for CuO as a function of tem¬ 
perature and vacancy concentration x for pure CuO and for 
CuO with X = 0.02 and x = 0.04. The horizontal dashed line 
indicates the cooling at constant entropy showing rather close 
quantitative agreement with the temperatures, indicated by 
vertical arrows with errorbars, measured within the spin dy¬ 
namics simulations. 



FIG. 9. Excitation intensity dependence of dynamics follow¬ 
ing on PD. The trajectory in time of S(q 2 ,t), averaged over 
128 replicas, and with a = 0.01. The initial temperature was 
T = 160 K, and the initial excitation intensity of vacancies 
was in the range [0.005, 0.040]. 


the lowering of T^i to a similar value requires a higher 
vacancy concentration, x si 0.03 [28]. This observation 
implies that the phase transition due to the lowering of 
Tjvi is possible when the system is kept in contact with 
the bath at the original temperature T = 160 K after 
the introduction of vacancies. To confirm this, we have 
carried out simulations at a = 0.01, shown by dashed 
lines in Fig. 7. In this case S'(qi,t) is further reduced 
after the fast initial drop, while S{q 2 ,t) increases. For 
vacancy concentration x = 0.02, the growth of S{q 2 ,t) is 
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still modest over the first ps, probably because the sys¬ 
tem is still to close to T;vi and the transition is broadened 
due to finite size effects. As shown in Fig. 9, a stronger 
response can be seen for doping up to a; = 0.04, for which 
also «S'(qi,t) decays completely and we are thus able to 
see the full transition on the time scale of our simulation. 
These results demonstrate that, opposed to the case of 
phonon-assisted MM excitation, the timescale at which 
the transition can occur is limited by the coupling to the 
bath and is therefore slower. 


C. Temperature induced dynamics 

For comparison, we also simulated direct heat-induced 
spin dynamics, mimicking a rapid rise of the phonon tem¬ 
perature by laser absorption. One can also take the bath 
to be excited e-h pairs in which case, in first approxima¬ 
tion, this effect should be added to the one of the PD 
computed above. The response of the spin system to 
an abrupt change of temperature at t = 0 is shown in 
Fig. 10. For simplicity, we consider a stepwise change of 
temperature which should give an indication of how fast 
the system react to heating. Unlike the case of phonon 
assisted MM excitation and PD, there is no instantaneous 
drop in the sublattice magnetization, so that S (qi, t) and 
S{q 2 ,t) are continuous at f = 0. This absence of instan¬ 
taneous decrease in the structure factor is in contrast 
with the diffraction data in [12]. We conclude that this 
mechanism is not playing an important role in those ex¬ 
periments. 

The speed of the response to heating is parametrically 
dependent on the strength of the damping parameter a 
[57]. Among the set of temperatures included in Fig. 10, 
the highest value of S(q 2 ,t) occur for T = 170 K. Unlike 
the case of phonon assisted MM excitation, the rate with 
which S'(q 2 , t) grows, is increasing for stronger pump flu- 
ence. 


VI. CONCLUSIONS 

In summary, we have explored how excitation with 
femtosecond pulses of near infrared laser light can drive 
a phase transition in CuO. On the basis of atomistic spin 
dynamics simulations we have found that photo-doping 
causes an ultrafast cooling of the spin system which 
can be understood as a magnetocaloric effect, where the 
collinear antiferromagnetic order drops rapidly and the 
spiral phase grows at a much slower rate due to the photo- 
induced reduction of Tati. Phonon-assisted MM excita¬ 
tion increase the temperature of the spin system on sub¬ 
picosecond timescales and a rapid drop of both collinear 
and spin spiral order is found, while only the spiral or¬ 
der grows gradually after this rapid drop. The phonon- 
assisted MM excitation is found to cause the fastest phase 
transition and we argue that this is related to the fact 
that only for this case the phase transition can proceed 


by pure exchange-driven spin dynamics [19, 55] alone. In 
addition, estimates of the quantum mechanical efficiency 
of inelastic light scattering was found to be too low to 
play any significant role in driving the phase transition, 
while a purely temperature driven dynamics cannot re¬ 
produce the initial abrupt reduction of the order param¬ 
eters observed in [12]. 

Finally, by comparing our results with experiments, we 
observe that only phonon-assisted MM absorption cap¬ 
tures the experimentally observed saturation with the 
pump fluence and the simultaneous initial drop of both 
collinear and spin-spiral order. On the other hand, in 
the experiments the timespan over which this simulta¬ 
neous drop takes place is a factor of 8 larger than in 
our simulations. There are several reasons for this dis¬ 
crepancy. First, the simulations assume an instantaneous 
excitation, whereas experimental laser pulses in [12] had 
a FWHM maximum of 40 fs. Moreover, in the experi¬ 
ment the laser pump pulse coincided with the band gap 
of the material, wherefore both PD decay and phonon 
assisted MM excitation are possible excitation mecha¬ 
nisms. On the level of ASD simulations, it is not straight¬ 
forward to realistically model of the interplay between 
these processes. In addition, we assume that concentra¬ 
tion of vacancies remains constant on the timescale con¬ 
sidered, while the eventual decay of the photo-excited 
holon-doublon pairs can also proceed through MM exci¬ 
tation [46-48]. Such decay can therefore drive the spin 
dynamics across the phase transition in a similar fashion 
as phonon-assisted MM excitation. Our simulations also 
show that these processes are not captured by a simple 
increase of the bath temperature. Similar to the absorp¬ 
tion of photons we expect that an initial state with MM 
excitation captures this effect more realistically. Clearly, 
modeling the simultaneous charge and spin dynamics 
goes beyond the semi-classical approach employed here 



FIG. 10. Dynamics following on a change of the heat bath 
temperature. The trajectory in time of S(qi,t) and S{q 2 ,t), 
averaged over 128 replicas, and with a = 0.01. The initial 
temperature was T = 160 K. 
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and would be very interesting to pursuit in future inves¬ 
tigations. It would be also interesting to repeat the ex¬ 
periments varying the laser wave-length to the iniddle-IR 
region to verify our predictions in the pure MM regime. 
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Appendix A: Measuring temperature 

In this appendix we compare two different methods to 
sample the temperature of a spin system: the method 
developed by Ma et al. [43] and our own construction 


based on energy histograms. The former has the advan¬ 
tage that a closed expression is used to measure the spin 
temperature (Eq. 16 of Ref. [43], in the following referred 
to as the MDSW equation), whereas an advantage of the 
latter is that the histograms also indicate how long it 
takes before the spin system is in internal equilibrium as 
signified by deviations from a Boltzmann distribution. 

In Figure 11(a) are shown in log-linear graphs the 
histograms for T = 100,150, 200, 250, and 300 K. For 
T = 150 K the bin count deviate from a straight line 
only for the highest energies, and the histogram can 
be fitted to a Boltzmann distribution with temperature 
Tfit = 150 ± 0.5 K. At T = 200 K the deviation is more 
pronounced and the fitted temperature is slightly under¬ 
estimated, Tfit ~ 196 K. At higher temperatures the de¬ 
viation is even more pronounced. 

For the bilinear in spin Hamiltonian + 

be. omitting the biquadratic interaction the 
MDSW equation measures the temperature accurately 
also in the paramagnetic phase, as shown in the lower 
inset of Fig. 11(b). However, for our full Hamiltonian 
the histogram method works better and does at a 
Tbath = 170 K undershoot with only 1 K. The deviations 
of the two methods from the heat bath temperatures are 
displayed in the upper inset of Fig. 11(b). For the simu¬ 
lations discussed in Sec. V, temperatures are in the range 
10 < Tfit < 170 K and the histogram method is expected 
to be sufficiently accurate. 
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FIG. 11. (a) Histograms (symbols) of energies Ei for spin 
systems in thermal equilibrium. The data for different tem¬ 
peratures have been shifted along the y-axis. The full lines 
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heat bath, dashed lines show the fitted temperatures, (b) 
The measured spin system temperature obtained with the 
histogram method (blue curve) or the MDSW formula (red 
curve). The upper inset displays the deviation of the mea¬ 
sure temperatures to the the bath temperature. The lower 
inset shows the measured spin system temperature for the 
bilinear in spin Hamiltonian J^mn- 
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